      subroutine bcphi(ibp1, jbp1, kbp1, iwid, jwid, kwid, wkl, wkr, wkf, wke, &
         phibc, phi, periodicx)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double
      use modify_com_M
      use boundary

!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02
!...Switches: -yf -x1
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer , intent(in) :: ibp1
      integer , intent(in) :: jbp1
      integer , intent(in) :: kbp1
      integer , intent(in) :: iwid
      integer , intent(in) :: jwid
      integer , intent(in) :: kwid
      real(double) , intent(in) :: wkl
      real(double) , intent(in) :: wkr
      real(double) , intent(in) :: wkf
      real(double) , intent(in) :: wke
      real(double) , intent(in) :: phibc
      logical , intent(in) :: periodicx
      real(double) , intent(inout) :: phi(*)
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: j, i, ijk, k, ijkl, ijkr, ijkb, ijkt, ijke, ijkf
      real(double) :: wkbloc,wktloc
!-----------------------------------------------
!
!     0=Dirichlet, 1-Neumann for wkf, wke, wkl, wkr, from the input file.
!
	if(periodicz) then
!
      do j = 1, jbp1+1
         do i = 1, ibp1+1
!
            ijke = 1 + (i - 1)*iwid + (j - 1)*jwid + kbp1*kwid
            ijkf = 1 + (i - 1)*iwid + (j - 1)*jwid
!
            phi(ijkf) = phi(ijke-kwid)
            phi(ijke) = phi(ijkf+kwid)
!
         end do
      end do

      else

!     k=kbp2 plane
!
      do j = 2, jbp1
         do i = 2, ibp1
!
            ijk = 1 + (i - 1)*iwid + (j - 1)*jwid + kbp1*kwid
!
!            phi(ijk) = (1. - wkf)*(2.*phibcf*phibc-phi(ijk-kwid)) + wkf*phi(ijk-kwid)
            phi(ijk) = (1. - wkf)*(phibcf*phibc) + wkf*phi(ijk-kwid)
!
         end do
      end do
!
!     k=1 plane (added by gl)
!
      do j = 2, jbp1
         do i = 2, ibp1
!
            ijk = 1 + (i - 1)*iwid + (j - 1)*jwid
!
!            phi(ijk) = (1. - wke)*(2.*phibce*phibc-phi(ijk+kwid)) + wke*phi(ijk+kwid)
            phi(ijk) = (1. - wke)*(phibce*phibc) + wke*phi(ijk+kwid)
!
         end do
      end do

      end if
!
!     commented out by gl cgl
!     minor axis of torus
!
!gl      do 4 j=2,jbp1
!
!      phibar=0.0
!      do 41 i=2,ibp1
!      ijk=1+(i-1)*iwid+(j-1)*jwid+kwid
!      phibar=phibar+phi(ijk)
!   41 continue
!
!      phibar=phibar/float(ibp1-1)
!
!
!dir$ ivdep
!gl      do 42 i=2,ibp1
!gl      ijk=1+(i-1)*iwid+(j-1)*jwid
!      phi(ijk)=2.*phibar-phi(ijk+kwid)
!gl      phi(ijk)=0.0
!gl   42 continue
!gl
!gl    4 continue
!
!      do 5 k=2,kbp1+1
      do k = 1, kbp1 + 1
!
!... i=1 & i=ibp2
         if (periodicx) then
!
!dir$ ivdep
            do j = 1, jbp1 + 1
               ijkl = 1 + (j - 1)*jwid + (k - 1)*kwid
               ijkr = 1 + ibp1*iwid + (j - 1)*jwid + (k - 1)*kwid
               phi(ijkr) = phi(ijkl+iwid)
               phi(ijkl) = phi(ijkr-iwid)
            end do
!
         else
!
!dir$ ivdep
            do j = 1, jbp1 + 1
               ijkl = 1 + (j - 1)*jwid + (k - 1)*kwid
               ijkr = 1 + ibp1*iwid + (j - 1)*jwid + (k - 1)*kwid
               phi(ijkr) = (1. - wkr)*(2.*phibc-phi(ijkr-iwid)) + wkr*phi(ijkr-iwid)
               phi(ijkl) = (1. - wkl)*(2.*phibc-phi(ijkl+iwid)) + wkl*phi(ijkl+iwid)
!      phi(ijkr)=2.*(1.-wkr)*phibc+(2.*wkr-1.)*phi(ijkr-iwid)
!      phi(ijkl)=2.*(1.-wkl)*phibc+(2.*wkl-1.)*phi(ijkl+iwid)
            end do
!
         endif
!
		if(periodicy) then
!dir$ ivdep
!... j=1 & j=jbp2,  periodic
         do i = 1, ibp1 + 1
            ijkb = 1 + (i - 1)*iwid + (k - 1)*kwid
            ijkt = 1 + (i - 1)*iwid + jbp1*jwid + (k - 1)*kwid
            phi(ijkb) = phi(ijkt-jwid)
            phi(ijkt) = phi(ijkb+jwid)
         end do
!
	    else
!... j=1 & j=jbp2,  periodic
if (bc_vtoctov) then
wkbloc=1d0
wktloc=1d0
else
wkbloc=wkb
wktloc=wkt
end if
         do i = 1, ibp1 + 1
            ijkb = 1 + (i - 1)*iwid + (k - 1)*kwid
            ijkt = 1 + (i - 1)*iwid + jbp1*jwid + (k - 1)*kwid
            phi(ijkb) = -(1. - wkbloc)*phi(ijkb+jwid) + wkbloc*phi(ijkb+jwid)
            phi(ijkt) = -(1. - wktloc)*phi(ijkt-jwid) + wktloc*phi(ijkt-jwid)
         end do
        end if
      end do
!
!
      return
      end subroutine bcphi



      subroutine bccz(nxp, nyp, nzp, fact, a)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double
      use modify_com_M

!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02
!...Switches: -yf -x1
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer , intent(in) :: nxp
      integer , intent(in) :: nyp
      integer , intent(in) :: nzp
      real(double) , intent(in) :: fact
      real(double) , intent(inout) :: a(nxp,nyp,*)
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: i, j
!-----------------------------------------------
!
      a(:,:,1) = fact*a(:,:,2)
      a(:,:,nzp) = fact*a(:,:,nzp-1)
!
      return
      end subroutine bccz



      subroutine bcvz(nxp, nyp, nzp, fact, a)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double
      use modify_com_M
      use boundary

!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02
!...Switches: -yf -x1
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer , intent(in) :: nxp
      integer , intent(in) :: nyp
      integer , intent(in) :: nzp
      real(double) , intent(in) :: fact
      real(double) , intent(inout) :: a(nxp,nyp,*)
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: i, j
!-----------------------------------------------
!
if (periodicz) then
      a(2:nxp,2:nyp,2) = a(2:nxp,2:nyp,2) + a(2:nxp,2:nyp,nzp)
      a(2:nxp,2:nyp,nzp) = a(2:nxp,2:nyp,2)
else
      a(2:nxp,2:nyp,2) = fact*a(2:nxp,2:nyp,2)
      a(2:nxp,2:nyp,nzp) = fact*a(2:nxp,2:nyp,nzp)
end if
!
      return
      end subroutine bcvz


!
      subroutine torusbc(nxp, nyp, nzp, cdlt, sdlt, strait, dz, x, y, z, &
         periodicx)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double
      use modify_com_M
      use boundary

!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02
!...Switches: -yf -x1
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer , intent(in) :: nxp
      integer , intent(in) :: nyp
      integer , intent(in) :: nzp
      real(double) , intent(in) :: cdlt
      real(double) , intent(in) :: sdlt
      real(double) , intent(in) :: strait
      real(double) , intent(in) :: dz
      logical :: periodicx
      real(double) , intent(inout) :: x(nxp,nyp,*)
      real(double) , intent(inout) :: y(nxp,nyp,*)
      real(double) , intent(inout) :: z(nxp,nyp,*)
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: k, i, j
!-----------------------------------------------
!
!
!     periodicity in the toroidal angle
!
    if(periodicy) then
!
      x(:,1,:nzp) = cdlt*x(:,nyp-1,:nzp) + sdlt*y(:,nyp-1,:nzp)
      y(:,1,:nzp) = (-sdlt*x(:,nyp-1,:nzp)) + cdlt*y(:,nyp-1,:nzp) - strait*dz
      z(:,1,:nzp) = z(:,nyp-1,:nzp)
!
      x(:,nyp,:nzp) = cdlt*x(:,2,:nzp) - sdlt*y(:,2,:nzp)
      y(:,nyp,:nzp) = sdlt*x(:,2,:nzp) + cdlt*y(:,2,:nzp) + strait*dz
      z(:,nyp,:nzp) = z(:,2,:nzp)

    else

      x(:,1,:nzp) = x(:,2,:nzp)
      y(:,1,:nzp) = y(:,2,:nzp)
      z(:,1,:nzp) = z(:,2,:nzp)
!
      x(:,nyp,:nzp) = x(:,nyp-1,:nzp)
      y(:,nyp,:nzp) = y(:,nyp-1,:nzp)
      z(:,nyp,:nzp) = z(:,nyp-1,:nzp)

 	end if
!
!
!
!     periodicity in the poloidal angle
!
      if (periodicx) then
!
!
         x(1,:,:nzp) = x(nxp-1,:,:nzp)
         y(1,:,:nzp) = y(nxp-1,:,:nzp)
         z(1,:,:nzp) = z(nxp-1,:,:nzp)
!
         x(nxp,:,:nzp) = x(2,:,:nzp)
         y(nxp,:,:nzp) = y(2,:,:nzp)
         z(nxp,:,:nzp) = z(2,:,:nzp)
!
!
      else
!
!
         x(1,:,:nzp) = x(2,:,:nzp)
         y(1,:,:nzp) = y(2,:,:nzp)
         z(1,:,:nzp) = z(2,:,:nzp)
!
         x(nxp,:,:nzp) = x(nxp-1,:,:nzp)
         y(nxp,:,:nzp) = y(nxp-1,:,:nzp)
         z(nxp,:,:nzp) = z(nxp-1,:,:nzp)
!
!
      endif
!
!
!     periodicity in the vertical direction
!
      if (periodicz) then
!
!
         x(:,:,1) = x(:,:,nzp-1)
         y(:,:,1) = y(:,:,nzp-1)
         z(:,:,1) = z(:,:,nzp-1)
!
         x(:,:,nzp) = x(:,:,2)
         y(:,:,nzp) = y(:,:,2)
         z(:,:,nzp) = z(:,:,2)
!
      else
!
         x(:,:,1) = x(:,:,2)
         y(:,:,1) = y(:,:,2)
         z(:,:,1) = z(:,:,2)
!
         x(:,:,nzp) = x(:,:,nzp-1)
         y(:,:,nzp) = y(:,:,nzp-1)
         z(:,:,nzp) = z(:,:,nzp-1)

!
      endif

!
      return
      end subroutine torusbc


      subroutine torusbg(nxp, nyp, nzp, istep, iota, cdlt, sdlt, strait, toroid&
         , rwall, rmaj, dz, x, y, z)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double
      use modify_com_M

!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02
!...Switches: -yf -x1
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer , intent(in) :: nxp
      integer , intent(in) :: nyp
      integer , intent(in) :: nzp
      real(double) , intent(in) :: cdlt
      real(double) , intent(in) :: sdlt
      real(double) , intent(in) :: strait
      real(double) , intent(in) :: toroid
      real(double) , intent(in) :: rwall
      real(double) , intent(in) :: rmaj
      real(double) , intent(in) :: dz
      integer , intent(inout) :: istep(*)
      real(double) , intent(in) :: iota(*)
      real(double) , intent(inout) :: x(nxp,nyp,*)
      real(double) , intent(inout) :: y(nxp,nyp,*)
      real(double) , intent(inout) :: z(nxp,nyp,*)
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: ibar, kbar, k, kt, i, i1, j
      real(double) :: delt, pi, dr, ws, wsr
!-----------------------------------------------
!
!     a routine to impose double periodicity
!     for toroidal geometry
!    N.B. Used only by grid generator
!
!
!     istep is recalculated as function of minor radius
!
      delt = acos(cdlt)
      pi = acos(-1.)
      ibar = nxp - 2
      kbar = nzp - 2
      dr = rwall/kbar
!
      do k = 2, nzp
         ws = toroid*sqrt(x(2,2,k)**2+y(2,2,k)**2) - rmaj + strait*x(2,2,k)
         wsr = sqrt(ws**2 + z(2,2,k)**2) + 1.E-20
         kt = ifix(wsr/dr) + 2
         istep(k) = ifix(iota(kt)*ibar+0.5)*delt/(2.*pi)
      end do
!
!
      do k = 2, nzp
!
!     periodicity in the toroidal angle
!
         do i = 1, nxp
!
            i1 = i + istep(k)
            if (i1 < 1) i1 = i1 + nxp - 2
            if (i1 > nxp) i1 = i1 - nxp + 2
!
            x(i1,1,k) = cdlt*x(i,nyp-1,k) + sdlt*y(i,nyp-1,k)
            y(i1,1,k) = (-sdlt*x(i,nyp-1,k)) + cdlt*y(i,nyp-1,k) - strait*dz
            z(i1,1,k) = z(i,nyp-1,k)
!
            x(i,nyp,k) = cdlt*x(i1,2,k) - sdlt*y(i1,2,k)
            y(i,nyp,k) = sdlt*x(i1,2,k) + cdlt*y(i1,2,k) + strait*dz
            z(i,nyp,k) = z(i1,2,k)
!
         end do
!
!
!     periodicity in the poloidal angle
!
!
         x(1,:,k) = x(nxp-1,:,k)
         y(1,:,k) = y(nxp-1,:,k)
         z(1,:,k) = z(nxp-1,:,k)
!
         x(nxp,:,k) = x(2,:,k)
         y(nxp,:,k) = y(2,:,k)
         z(nxp,:,k) = z(2,:,k)
!
!
      end do
!
      return
      end subroutine torusbg



      subroutine bcperv(nxp, nyp, nzp, iwid, jwid, kwid, sdlt, cdlt, ax, ay, az&
         , periodicx)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double
      use modify_com_M
      use boundary

!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02
!...Switches: -yf -x1
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer , intent(in) :: nxp
      integer , intent(in) :: nyp
      integer , intent(in) :: nzp
      integer , intent(in) :: iwid
      integer , intent(in) :: jwid
      integer , intent(in) :: kwid
      real(double) , intent(in) :: sdlt
      real(double) , intent(in) :: cdlt
      logical :: periodicx
      real(double) , intent(inout) :: ax(*)
      real(double) , intent(inout) :: ay(*)
      real(double) , intent(inout) :: az(*)
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: k, j, ijk1, ijknxp, i, ijknyp, ijknzp
!-----------------------------------------------
!
!     set reference vectors for i=1 and i=nxp
!
      do k = 1, nzp
!
         if (periodicx) then
!
            do j = 1, nyp
!
               ijk1 = 1 + (j - 1)*jwid + (k - 1)*kwid
               ijknxp = 1 + (nxp - 1)*iwid + (j - 1)*jwid + (k - 1)*kwid
!
               ax(ijk1) = ax(ijknxp-iwid)
               ax(ijknxp) = ax(ijk1+iwid)
!
               ay(ijk1) = ay(ijknxp-iwid)
               ay(ijknxp) = ay(ijk1+iwid)
!
               az(ijk1) = az(ijknxp-iwid)
               az(ijknxp) = az(ijk1+iwid)
!
            end do
!
         else
!
            do j = 1, nyp
!
               ijk1 = 1 + (j - 1)*jwid + (k - 1)*kwid
               ijknxp = 1 + (nxp - 1)*iwid + (j - 1)*jwid + (k - 1)*kwid
!
               ax(ijk1) = ax(ijk1+iwid)
!      ax(ijknxp)=ax(ijknxp-iwid)
!
               ay(ijk1) = ay(ijk1+iwid)
!      ay(ijknxp)=ay(ijknxp-iwid)
!
               az(ijk1) = az(ijk1+iwid)
!      az(ijknxp)=az(ijknxp-iwid)
!
            end do
!
         endif
!
!     set reference vectors for j=1 and j=nyp
!
        if(periodicy) then
!
         do i = 1, nxp
!
            ijk1 = 1 + (i - 1)*iwid + (k - 1)*kwid
            ijknyp = 1 + (i - 1)*iwid + (nyp - 1)*jwid + (k - 1)*kwid
!
            ax(ijk1) = ax(ijknyp-jwid)*cdlt + ay(ijknyp-jwid)*sdlt
            ay(ijk1) = (-ax(ijknyp-jwid)*sdlt) + ay(ijknyp-jwid)*cdlt
            az(ijk1) = az(ijknyp-jwid)
!
            ax(ijknyp) = ax(ijk1+jwid)*cdlt - ay(ijk1+jwid)*sdlt
            ay(ijknyp) = ax(ijk1+jwid)*sdlt + ay(ijk1+jwid)*cdlt
            az(ijknyp) = az(ijk1+jwid)
!
         end do
!
	   else
!
         do i = 1, nxp
!
            ijk1 = 1 + (i - 1)*iwid + (k - 1)*kwid
!
            ax(ijk1) = ax(ijk1+jwid)
            ay(ijk1) = ay(ijk1+jwid)
            az(ijk1) = az(ijk1+jwid)
!
         end do
!
	   end if
!
      end do
!
!     Set reference vectors for k=1 and k=nzp
!
      if (periodicz) then
           do i = 1, nxp
!
         do j = 1, nyp
!
            ijk1 = 1 + (i - 1)*iwid + (j - 1)*jwid
            ijknzp = 1 + (i - 1)*iwid + (j - 1)*jwid + (nzp - 1)*kwid
!
            ax(ijk1) = ax(ijknzp-kwid)
            ax(ijknzp) = ax(ijk1+kwid)
            ay(ijk1) = ay(ijknzp-kwid)
            ay(ijknzp) = ay(ijk1+kwid)
            az(ijk1) = az(ijknzp-kwid)
            az(ijknzp) = az(ijk1+kwid)
!
         end do
      end do

      else

      do i = 1, nxp
!
         do j = 1, nyp
!
            ijk1 = 1 + (i - 1)*iwid + (j - 1)*jwid
            ijknzp = 1 + (i - 1)*iwid + (j - 1)*jwid + (nzp - 1)*kwid
!
            ax(ijk1) = ax(ijk1+kwid)
!      ax(ijknzp)=ax(ijknzp-kwid)
!      ax(ijk1)=0.
!      ax(ijknzp)=0.
!
            ay(ijk1) = ay(ijk1+kwid)
!      ay(ijknzp)=ay(ijknzp-kwid)
!      ay(ijk1)=0.
!      ay(ijknzp)=0.
!
            az(ijk1) = az(ijk1+kwid)
!      az(ijknzp)=az(ijknzp-kwid)
!      az(ijk1)=0.
!      az(ijknzp)=0.
!
         end do
      end do
      end if

      return
      end subroutine bcperv



      subroutine torusbcv(nxp, nyp, nzp, cdlt, sdlt, strait, dz, x, y, z, &
         periodicx)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double
      use modify_com_M
      use boundary

!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02
!...Switches: -yf -x1
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer , intent(in) :: nxp
      integer , intent(in) :: nyp
      integer , intent(in) :: nzp
      real(double) , intent(in) :: cdlt
      real(double) , intent(in) :: sdlt
      real(double) , intent(in) :: strait
      real(double) , intent(in) :: dz
      logical , intent(in) :: periodicx
      real(double) , intent(inout) :: x(nxp,nyp,*)
      real(double) , intent(inout) :: y(nxp,nyp,*)
      real(double) , intent(inout) :: z(nxp,nyp,*)
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: k, i, j
!-----------------------------------------------
!
!
!
!     periodicity in the toroidal angle
!
!

	if(periodicy) then

      x(2:nxp,2,2:nzp) = cdlt*x(2:nxp,nyp,2:nzp) + sdlt*y(2:nxp,nyp,2:nzp) + x(&
         2:nxp,2,2:nzp)
      y(2:nxp,2,2:nzp) = (-sdlt*x(2:nxp,nyp,2:nzp)) + cdlt*y(2:nxp,nyp,2:nzp)&
          + y(2:nxp,2,2:nzp) - strait*dz
      z(2:nxp,2,2:nzp) = z(2:nxp,nyp,2:nzp) + z(2:nxp,2,2:nzp)
!
      x(2:nxp,nyp,2:nzp) = cdlt*x(2:nxp,2,2:nzp) - sdlt*y(2:nxp,2,2:nzp)
      y(2:nxp,nyp,2:nzp) = sdlt*x(2:nxp,2,2:nzp) + cdlt*y(2:nxp,2,2:nzp) + &
         strait*dz
      z(2:nxp,nyp,2:nzp) = z(2:nxp,2,2:nzp)

	else

      x(2:nxp,2,2:nzp) = 2d0 * x(2:nxp,2,2:nzp)
      y(2:nxp,2,2:nzp) = 2d0 * y(2:nxp,2,2:nzp)
      z(2:nxp,2,2:nzp) = 2d0 * z(2:nxp,2,2:nzp)
!
      x(2:nxp,nyp,2:nzp) = 2d0 * x(2:nxp,nyp,2:nzp)
      y(2:nxp,nyp,2:nzp) = 2d0 * y(2:nxp,nyp,2:nzp)
      z(2:nxp,nyp,2:nzp) = 2d0 * z(2:nxp,nyp,2:nzp)

	end if
!
!
!
!     periodicity in the poloidal angle
!
      if (periodicx) then
!
!
         x(2,2:nyp,2:nzp) = x(nxp,2:nyp,2:nzp) + x(2,2:nyp,2:nzp)
         y(2,2:nyp,2:nzp) = y(nxp,2:nyp,2:nzp) + y(2,2:nyp,2:nzp)
         z(2,2:nyp,2:nzp) = z(nxp,2:nyp,2:nzp) + z(2,2:nyp,2:nzp)
!
         x(nxp,2:nyp,2:nzp) = x(2,2:nyp,2:nzp)
         y(nxp,2:nyp,2:nzp) = y(2,2:nyp,2:nzp)
         z(nxp,2:nyp,2:nzp) = z(2,2:nyp,2:nzp)
!
!
      else
!
!
         x(2,2:nyp,2:nzp) = 2.*x(2,2:nyp,2:nzp)
         y(2,2:nyp,2:nzp) = 2.*y(2,2:nyp,2:nzp)
         z(2,2:nyp,2:nzp) = 2.*z(2,2:nyp,2:nzp)
!
         x(nxp,2:nyp,2:nzp) = 2.*x(nxp,2:nyp,2:nzp)
         y(nxp,2:nyp,2:nzp) = 2.*y(nxp,2:nyp,2:nzp)
         z(nxp,2:nyp,2:nzp) = 2.*z(nxp,2:nyp,2:nzp)
!
!
      endif
!
!
!
      return
      end subroutine torusbcv




      subroutine torusbc_scal(i1, i2, j1, j2, k1, k2, iwid, jwid, kwid, pars, &
         qdnv, periodicx)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double
      use modify_com_M
      use boundary

!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02
!...Switches: -yf -x1
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer , intent(in) :: i1
      integer , intent(in) :: i2
      integer , intent(in) :: j1
      integer , intent(in) :: j2
      integer , intent(in) :: k1
      integer , intent(in) :: k2
      integer , intent(in) :: iwid
      integer , intent(in) :: jwid
      integer , intent(in) :: kwid
      real(double) , intent(in) :: pars
      logical , intent(in) :: periodicx
      real(double) , intent(inout) :: qdnv(*)
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: k, i, ijkb, ijkt, j, ijkl, ijkr, ijk
      real(double) :: qdvsum, fnorm, rnorm
!-----------------------------------------------
!
      do k = k1, k2

       if(periodicy) then

         do i = i1, i2
            ijkb = 1 + (i - 1)*iwid + (j1 - 1)*jwid + (k - 1)*kwid
            ijkt = 1 + (i - 1)*iwid + (j2 - 1)*jwid + (k - 1)*kwid
!
            qdnv(ijkt) = pars*qdnv(ijkt) + qdnv(ijkb)
            qdnv(ijkb) = qdnv(ijkt)
!      qdnv(ijkb-jwid)=qdnv(ijkt-jwid)
!
         end do
!
       else

         do i = i1, i2
            ijkb = 1 + (i - 1)*iwid + (j1 - 1)*jwid + (k - 1)*kwid
            ijkt = 1 + (i - 1)*iwid + (j2 - 1)*jwid + (k - 1)*kwid
!
            qdnv(ijkt) = pars*2d0* qdnv(ijkt)
            qdnv(ijkb) = pars*2d0* qdnv(ijkb)
!
         end do
!
       end if
!
         if (periodicx) then
!
            do j = j1, j2
               ijkl = 1 + (i1 - 1)*iwid + (j - 1)*jwid + (k - 1)*kwid
               ijkr = 1 + (i2 - 1)*iwid + (j - 1)*jwid + (k - 1)*kwid
!
               qdnv(ijkr) = pars*qdnv(ijkr) + qdnv(ijkl)
               qdnv(ijkl) = qdnv(ijkr)
!      qdnv(ijkl-iwid)=qdnv(ijkr-iwid)
!
            end do
!
         else
!
            do j = j1, j2
               ijkl = 1 + (i1 - 1)*iwid + (j - 1)*jwid + (k - 1)*kwid
               ijkr = 1 + (i2 - 1)*iwid + (j - 1)*jwid + (k - 1)*kwid
!
               qdnv(ijkr) = pars*2.*qdnv(ijkr)
               qdnv(ijkl) = pars*2.*qdnv(ijkl)
!      qdnv(ijkr)=pars*qdnv(ijkr)
!      qdnv(ijkl)=pars*qdnv(ijkl)
!
            end do
!
         endif
!
      end do
!     WARNING THIS SUBROUTINE ONLY WORKS FOR CARTESIAN GEOMETRY
      return
!
      end subroutine torusbc_scal

      subroutine boundary_b_centered(nxp,nyp,nzp, dx,dy,dz, dt, &
         bxn, byn, bzn, periodicx, neumann)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double
      use modify_com_M
      use boundary
!
! Routine to impose bopundary conditions on the magnetic field B
! defined on cell centres.
!
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer  :: nxp,nyp,nzp
      logical  :: periodicx
      real(double)  :: dx,dy,dz,dt, neumann
      real(double)  :: bxn(nxp, nyp, nzp)
      real(double)  :: byn(nxp, nyp, nzp)
      real(double)  :: bzn(nxp, nyp, nzp)
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: i,j,k
!-----------------------------------------------
!-----------------------------------------------
!
!
!	Boundaries in Y
!
    if (periodicy) then
	do i=1,nxp
	do k=1,nzp
	   bxn(i,1,k)=bxn(i,nyp-1,k)
	   byn(i,1,k)=byn(i,nyp-1,k)
	   bzn(i,1,k)=bzn(i,nyp-1,k)
	   bxn(i,nyp,k)=bxn(i,2,k)
	   byn(i,nyp,k)=byn(i,2,k)
	   bzn(i,nyp,k)=bzn(i,2,k)
	enddo
	enddo

	else

    if(bcMy.eq.1) then

	!apply conductor
	do k=1,nzp
	do i=1,nxp
	   bxn(i,1,k)=bxn(i,2,k)
	   byn(i,1,k)= - byn(i,2,k)
	   bzn(i,1,k)= bzn(i,2,k)
	   bxn(i,nyp,k)=bxn(i,nyp-1,k)
	   byn(i,nyp,k)= - byn(i,nyp-1,k)
	   bzn(i,nyp,k)=bzn(i,nyp-1,k)
	enddo
	enddo

	elseif(bcMy.eq.2) then
	!apply mirror
	do k=1,nzp
	do i=1,nxp
	   bxn(i,1,k)=0d0
	   byn(i,1,k)= byn(i,2,k)
	   bzn(i,1,k)= 0d0
	   bxn(i,nyp,k)=0d0
	   byn(i,nyp,k)= byn(i,nyp-1,k)
	   bzn(i,nyp,k)=0d0
	enddo
	enddo


    elseif(bcMy.eq.0) then

        !apply open

	do k=1,nzp
	do i=1,nxp
	   bxn(i,1,k)=neumann*bxn(i,2,k) - (1d0-neumann)*bxn(i,2,k)
	   byn(i,1,k)= byn(i,2,k)
	   bzn(i,1,k)= neumann*bzn(i,2,k) - (1d0-neumann)*bzn(i,2,k)
	   bxn(i,nyp,k)=neumann*bxn(i,nyp-1,k) - (1d0-neumann)*bxn(i,nyp-1,k)
	   byn(i,nyp,k)= byn(i,nyp-1,k)
	   bzn(i,nyp,k)=neumann*bzn(i,nyp-1,k) - (1d0-neumann)*bzn(i,nyp-1,k)
	enddo
	enddo

    else

        !apply open ... Pritchett ...
	do k=1,nzp
	do i=1,nxp
	   bxn(i,1,k)=neumann*bxn(i,2,k) - (1d0-neumann)*bxn(i,2,k)
	   byn(i,1,k)= byn(i,2,k)
	   bzn(i,1,k)= neumann*bzn(i,2,k)! - (1d0-neumann)*bzn(i,2,k)
	   bxn(i,nyp,k)=neumann*bxn(i,nyp-1,k) - (1d0-neumann)*bxn(i,nyp-1,k)
	   byn(i,nyp,k)= byn(i,nyp-1,k)
	   bzn(i,nyp,k)=neumann*bzn(i,nyp-1,k)! - (1d0-neumann)*bzn(i,nyp-1,k)
	enddo
	enddo

	end if

	end if
!
!	Boundaries in X
!
	if(periodicx) then

	do j=1,nyp
	do k=1,nzp
	   bxn(1,j,k)=bxn(nxp-1,j,k)
	   byn(1,j,k)=byn(nxp-1,j,k)
	   bzn(1,j,k)=bzn(nxp-1,j,k)
	   bxn(nxp,j,k)=bxn(2,j,k)
	   byn(nxp,j,k)=byn(2,j,k)
	   bzn(nxp,j,k)=bzn(2,j,k)
	enddo
	enddo

	else

	if(bcMx.eq.1) then

	!apply conductor
	do j=1,nyp
	do k=1,nzp
	   bxn(1,j,k)=-bxn(2,j,k)
	   byn(1,j,k)=byn(2,j,k)
	   bzn(1,j,k)=bzn(2,j,k)
	   bxn(nxp,j,k)=-bxn(nxp-1,j,k)
	   byn(nxp,j,k)=byn(nxp-1,j,k)
	   bzn(nxp,j,k)=bzn(nxp-1,j,k)
	enddo
	enddo

	elseif(bcMx.eq.2) then
!        if(bcMx.eq.0.and.neumann>0.5) goto 123 ! open and for B not curlE
	!apply mirror
	do j=1,nyp
	do k=1,nzp
	   bxn(1,j,k)=bxn(2,j,k)
	   byn(1,j,k)=0d0
	   bzn(1,j,k)=0d0
	   bxn(nxp,j,k)=bxn(nxp-1,j,k)
	   byn(nxp,j,k)=0d0
	   bzn(nxp,j,k)=0d0
	enddo
	enddo

    elseif(bcMx.eq.0) then
    !apply open
    do j=1,nyp
	do k=1,nzp
	   bxn(1,j,k)=bxn(2,j,k)
	   byn(1,j,k)=neumann*byn(2,j,k) - (1d0-neumann)*byn(2,j,k)
	   bzn(1,j,k)=neumann*bzn(2,j,k) - (1d0-neumann)*bzn(2,j,k)
	   bxn(nxp,j,k)=bxn(nxp-1,j,k)
	   byn(nxp,j,k)=neumann*byn(nxp-1,j,k) - (1d0-neumann)*byn(nxp-1,j,k)
	   bzn(nxp,j,k)=neumann*bzn(nxp-1,j,k) - (1d0-neumann)*bzn(nxp-1,j,k)
	enddo
	enddo

        elseif(bcMx.eq.-1) then
        !apply open
        do j=1,nyp
        do k=1,nzp
           bxn(1,j,k)=bxn(2,j,k)
           byn(1,j,k)=byn(2,j,k)
           bzn(1,j,k)= -bzn(2,j,k)
           bxn(nxp,j,k)=bxn(nxp-1,j,k)
           byn(nxp,j,k)=byn(nxp-1,j,k)
           bzn(nxp,j,k)= -bzn(nxp-1,j,k)
        enddo
        enddo

        elseif(bcMx.eq.-2) then
        !apply open
        do j=1,nyp
        do k=1,nzp
           bxn(1,j,k)=(5d0*bxn(2,j,k)-4d0*bxn(3,j,k)+bxn(4,j,k))/2d0
           byn(1,j,k)=byn(2,j,k)
           bzn(1,j,k)=bzn(2,j,k)
           bxn(nxp,j,k)=(5d0*bxn(nxp-1,j,k)-4d0*bxn(nxp-2,j,k)+bxn(nxp-3,j,k))/2d0
           byn(nxp,j,k)=byn(nxp-1,j,k)
           bzn(nxp,j,k)=bzn(nxp-1,j,k)
        enddo
        enddo

        else
        !apply open, Pritchett, actually not a real open BC..
        do j=1,nyp
	do k=1,nzp
	   bxn(1,j,k)=bxn(2,j,k)
	   byn(1,j,k)=neumann*byn(2,j,k)! - (1d0-neumann)*byn(2,j,k)
	   bzn(1,j,k)=neumann*bzn(2,j,k) - (1d0-neumann)*bzn(2,j,k)
	   bxn(nxp,j,k)=bxn(nxp-1,j,k)
	   byn(nxp,j,k)=neumann*byn(nxp-1,j,k)! - (1d0-neumann)*byn(nxp-1,j,k)
	   bzn(nxp,j,k)=neumann*bzn(nxp-1,j,k) - (1d0-neumann)*bzn(nxp-1,j,k)
	enddo
	enddo
	end if

	end if
!
!	Boundaries in Z
!
123	if(periodicz) then

	do j=1,nyp
	do i=1,nxp
	   bxn(i,j,1)=bxn(i,j,nzp-1)
	   byn(i,j,1)=byn(i,j,nzp-1)
	   bzn(i,j,1)=bzn(i,j,nzp-1)
	   bxn(i,j,nzp)=bxn(i,j,2)
	   byn(i,j,nzp)=byn(i,j,2)
	   bzn(i,j,nzp)=bzn(i,j,2)
	enddo
	enddo

	else

	if(bcMz.eq.1) then

	!apply conductor
	do j=1,nyp
	do i=1,nxp
	   bxn(i,j,1)=bxn(i,j,2)
	   byn(i,j,1)=byn(i,j,2)
	   bzn(i,j,1)= - bzn(i,j,2) ! for newton challenge this ought to be +? No.
	   bxn(i,j,nzp)=bxn(i,j,nzp-1)
	   byn(i,j,nzp)=byn(i,j,nzp-1)
	   bzn(i,j,nzp)= - bzn(i,j,nzp-1) ! for newton challenge this ought to be +? No.
	enddo
	enddo

	elseif(bcMz.eq.2) then
!	if(bcMz.eq.0.and.neumann>0.5) return ! open and for B not curlE
	!apply mirror
	do j=1,nyp
	do i=1,nxp
	   bxn(i,j,1)=0d0
	   byn(i,j,1)=0d0
	   bzn(i,j,1)=bzn(i,j,2)
	   bxn(i,j,nzp)=0d0
	   byn(i,j,nzp)=0d0
	   bzn(i,j,nzp)=bzn(i,j,nzp-1)
	enddo
	enddo

    elseif(bcMz.eq.0) then

        !apply open
	do j=1,nyp
	do i=1,nxp
	   bxn(i,j,1)=neumann*bxn(i,j,2) - (1d0-neumann)*bxn(i,j,2)
	   byn(i,j,1)=neumann*byn(i,j,2) - (1d0-neumann)*byn(i,j,2)
	   bzn(i,j,1)=bzn(i,j,2)
	   bxn(i,j,nzp)=neumann*bxn(i,j,nzp-1) - (1d0-neumann)*bxn(i,j,nzp-1)
	   byn(i,j,nzp)=neumann*byn(i,j,nzp-1) - (1d0-neumann)*byn(i,j,nzp-1)
	   bzn(i,j,nzp)=bzn(i,j,nzp-1)
	enddo
	enddo

    else

        !apply open ... Pritchett ...
	do j=1,nyp
	do i=1,nxp
	   bxn(i,j,1)=neumann*bxn(i,j,2) - (1d0-neumann)*bxn(i,j,2)
	   byn(i,j,1)=neumann*byn(i,j,2)! - (1d0-neumann)*byn(i,j,2)
	   bzn(i,j,1)=bzn(i,j,2)
	   bxn(i,j,nzp)=neumann*bxn(i,j,nzp-1) - (1d0-neumann)*bxn(i,j,nzp-1)
	   byn(i,j,nzp)=neumann*byn(i,j,nzp-1)! - (1d0-neumann)*byn(i,j,nzp-1)
	   bzn(i,j,nzp)=bzn(i,j,nzp-1)
	enddo
	enddo

	end if

	end if
	end subroutine boundary_b_centered

     subroutine boundary_b_vertex(nxp,nyp,nzp, dx,dy,dz, dt, &
         bxn, byn, bzn, periodicx, neumann)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double
      use modify_com_M
      use boundary
!
! Routine to impose bopundary conditions on the magnetic field B
! defined on cell centres.
!
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer  :: nxp,nyp,nzp
      logical  :: periodicx
      real(double)  :: dx,dy,dz,dt, neumann
      real(double)  :: bxn(nxp, nyp, nzp)
      real(double)  :: byn(nxp, nyp, nzp)
      real(double)  :: bzn(nxp, nyp, nzp)
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: i,j,k
      real(double) :: density,temperature,flvx,flvy,flvz,bfldx,bfldy,bfldz
!-----------------------------------------------
!-----------------------------------------------

    if(bcMz.eq.1.and..not.periodicz) then

    !apply conductor
    do j=1,nyp
    do i=1,nxp
       k=2
       call boundary_state(i,j,k,density,temperature,flvx,flvy,flvz,bfldx,bfldy,bfldz)
       !if(i==1.and.j==1) write(*,*)'boundary B', bfldx,bfldy,bfldz,density,temperature
       bxn(i,j,2)=bfldx
       byn(i,j,2)=bfldy
       bzn(i,j,2)=bfldz
       bxn(i,j,1)=bfldx
       byn(i,j,1)=bfldy
       bzn(i,j,1)=bfldz

       k=nzp-1
       call boundary_state(i,j,k,density,temperature,flvx,flvy,flvz,bfldx,bfldy,bfldz)
       !if(i==1.and.j==1) write(*,*)'boundary B', bfldx,bfldy,bfldz,density,temperature
       bxn(i,j,nzp)=bfldx
       byn(i,j,nzp)=bfldy
       bzn(i,j,nzp)=bfldz
    enddo
    enddo

    end if

    end subroutine boundary_b_vertex
